Human pancreatic islet miRNA-mRNA networks of altered miRNAs due to glycemic status

Summary MicroRNAs (miRNAs) are short non-coding RNAs that regulate gene expression via mRNA targeting, playing important roles in the pancreatic islets. We aimed to identify molecular pathways and genomic regulatory regions associated with altered miRNA expression due to glycemic status, which could contribute to the development of type 2 diabetes (T2D). To this end, miRNAs were identified by a combination of differential miRNA expression and correlation analysis in human islet samples from donors with normal and elevated blood glucose levels. Analysis and clustering of highly correlated, experimentally validated gene targets of these miRNAs revealed two islet-specific clusters, which were associated with key aspects of islet functions and included a high number of T2D-related genes. Finally, cis-eQTLs and public GWAS data integration uncovered suggestive genomic signals of association with insulin secretion and T2D. The miRNA-driven network-based approach presented in this study contributes to a better understanding of impaired insulin secretion in T2D pathogenesis.


INTRODUCTION
miRNAs are short (z19-23nt) endogenous non-coding RNAs which most of the time are responsible for silencing genes through inhibition of their translation or destabilization of their target mRNA molecule (Bartel, 2009). However, some cases of upregulation of the miRNA target genes have also been reported (Vasudevan et al., 2007). They have important regulatory roles in various biological processes, from specifying cell identity during development, to fine-tuning cellular functions in response to environmental stimuli. Indeed, it is estimated that 60% of mammalian protein-coding genes are posttranscriptionally regulated by miRNAs (Friedman et al., 2009). Unsurprisingly, perturbed miRNA expression in different tissues has become a recognized feature in human disease pathophysiology (Paul et al., 2018).
Type 2 diabetes (T2D) is characterized by impairment in the regulation and utilization of blood glucose as an energy source in the body. Contributing factors of the pathogenesis is the combination of insulin resistance, in which target tissues lose the capacity of properly responding to insulin, and the dysfunction of pancreatic b-cells, which secrete suboptimal levels of insulin (DeFronzo et al., 2015;Lebovitz, 1999). The functional impairment of glucagon-secreting pancreatic a-cells could also be responsible for the progression of the disease, because glucagon counteracts the effects of insulin and is crucial for maintaining the glucose homeostasis (Wendt and Eliasson, 2020). The importance of miRNAs in the pancreatic islet cells has been unequivocally shown in knockout mouse models of specific miRNAs or via global ablation of miRNA expression by deleting the Dicer1 enzyme that is involved in miRNA maturation (Kalis et al., 2011;Lynn et al., 2007;Martinez-Sanchez et al., 2015;Melkman-Zehavi et al., 2011). Such studies demonstrate that miRNAs can be implicated in events leading to T2D onset. In T2D, b-cell compensation is important, with some miRNAs exerting compensatory effects and others impacting insulin secretion through miRNA-mediated dysfunction (Eliasson and Regazzi, 2020). Among specific islet miRNAs shown to be implicated in dysfunctional b-cell function are miR-7, miR130a/130b, miR-152, miR-184, and the miR-200 family (Eliasson and Regazzi, 2020). In islets, most studies so far have investigated the role of single miRNAs and have not used a network approach of differentially expressed miRNAs and their targets. However, we and others suggest that islet function is regulated by miRNA groups rather than single miRNAs (Eliasson and Esguerra, 2020;Wong et al., 2021).
Genome-wide association studies (GWAS) have been effective in connecting genomic regions to various phenotypes, including T2D and related glycemic traits (Mahajan et al., 2014;Morris et al., 2012). The fact that individual T2D-associated variants have a small contribution to the overall disease risk  and that the vast majority of them are located in non-coding regions (Mahajan et al., 2018), makes the association between the variants with specific proteins and pathways complicated. Nevertheless, the identification of disease-specific GWAS variants that are related to the expression levels of proximal genes, or expression quantitative trait loci (eQTL), enables formulation of hypotheses about specific proteins/pathways (Westra and Franke, 2014).
In this study, we aimed to discover human islet miRNAs that are affected due to changes in the glycemic status and investigate their targeting properties by retrieving a list of validated gene targets that are expressed in human islets. Moreover, we identified pathways of the differentially expressed miRNAs and their targets, their functional annotation, and their correlation with different phenotypic traits. Finally, we investigated genomic regions that are associated with miRNA expression via cis-eQTL analysis and found SNPs linked to T2D and insulin secretion traits.

Global miRNA profiling of human pancreatic islets
Human islets from 18 donors (Table S1) were profiled for 840 miRNAs using Exiqon's LNA (locked nucleic acid)-based array platform. In total, 470 miRNAs were expressed in our islet samples (Table S2).
The donors were divided into three groups according to their HbA1c levels. The HbA1c values of the 18 donors ranges from 4.6% to 7%, seven had normal glucose tolerance (NGT) (HbA1c<6%), six had impaired glucose tolerance (IGT) (HBA1c range: 6.0%-6.4%), while five of them were also diagnosed with T2D (HbA1c range: 6.2%-7.0%) (Table S1). Glucose-stimulated insulin secretion (GSIS) data were available for four NGT, five IGT, and four T2D of our islet preparations. Data indicate a decreased trend regarding the function of the IGT and T2D islets compared to the NGT islets ( Figure 1A).
By analyzing probes that were captured by all donors and excluding those that corresponded to more than one human miRNAs, the remaining 269 miRNAs showed a large variation regarding their expression levels with only a few miRNAs being highly expressed ( Figure S1). Among the 15 most highly expressed miRNAs, we observed both miRNAs that were previously described as highly abundant in human islets or b-cells such as miR-375, miR-141-3p, and members of miR-29, miR-200, and miR-7 families (Kameswaran et al., 2014;Locke et al., 2014) and miRNAs not previously shown to be highly abundant, such as miR-1246, miR-1290, and miR-1908 ( Figure 1B). Despite the slight differences in their expression levels and ranking, 11 out of the 15 most highly expressed miRNAs across all donors were still abundant in islets derived from NGT, IGT, and T2D donors ( Figure 1B).
Using transcriptome data from the same donors derived from a previous study , we sought to reveal correlations between miRNA expression and gene expression of the major islet hormones, namely insulin, glucagon, and somatostatin. The analyses showed 19 miRNAs to be significantly correlated with insulin, 7 with glucagon, and 24 with somatostatin expression (p<0.05), while miR-27b-5p was correlated with the expression of both insulin and glucagon and miRNAs miR-21-3p and miR-497-5p to both glucagon and somatostatin (Table S3).

Altered expression of miRNAs due to glycemic status
The workflow for the identification of islet miRNAs with differential expression due to glycemic status of the donors can be seen in Figure 2. First, we matched the NGT donors with the IGT/T2D donors for age, gender, and BMI. The classification enabled the differential expression analysis of miRNAs between the two paired groups using Significance Analysis of Microarrays (SAM) based on false discovery rate (FDR) statistics within the TM4 platform (Saeed et al., 2003) (Figure 2, Table S4). The analysis revealed 37 upregulated and 26 downregulated miRNAs in islets from the IGT/T2D donors compared with NGT donors (FDR<0.1). Among the upregulated miRNAs were miR-1275, miR-32-3p, and miR-130b-5p/3p and among the downregulated we found miR-126-3p, miR-7-5p, and miR-200a.
Next, individual miRNA levels in all 18 donors were associated with HbA1c levels to assess how the in vivo long-term glucose exposure affects islet miRNA regulation. For this, we used a linear regression model and  (Table S4). Specifically, we could detect miR-23b, miR-7, and miR-484 to be correlated with HbA1c. All three are miRNAs previously shown to be regulated by glucose (Tang et al., 2009).

Gene target analysis of miRNAs with altered expression due to glycemic status
We were curious to explore the pathways in which these two miRNA datasets (Up-Pos and Down-Neg) are involved and, consequently, altered in response to glycemic status. To this end, we first determined the iScience Article gene targets of each miRNA by obtaining a list of experimentally validated gene targets for our miRNAs of interest. Since target validation of large miRNA sets has been facilitated by the introduction of highthroughput methods, validated gene targets of each miRNA were derived from Tarbase v.8 (Karagkouni et al., 2018) and miRTarBase v.7 (Chou et al., 2018) databases (Table S5). The nine Up-Pos miRNAs were associated with 4,046 unique genes forming 5,232 miRNA-mRNA interactions and the 22 Down-Neg miR-NAs were associated with 12,289 unique genes forming 41,950 miRNA-mRNA interactions. The higher number of interactions can be explained by the fact that different miRNAs can target the same gene and a single gene can be targeted by multiple miRNAs (Eliasson and Esguerra, 2020). We could show that the same number of random miRNAs could not achieve similar target/interaction numbers (Figure S2A). The Down-Neg miRNAs shared 29% of their gene targets with the Up-Pos miRNAs ( Figure 3A).
Considering that the validated gene target lists are collected from experiments in different cell types or tissues, we decided to retain only genes that are expressed in human islets for further investigation. To this end, we made use of the published human islet transcriptome data of 188 donors . Expressed genes were considered those with at least three normalized counts in at least 80% of the samples, leaving a total of 11,689 out of the 60,517 identified genes (Table S6). We showed that a high percentage of the validated gene targets of each miRNA are expressed in the human islets (miRNA with lowest overlap 68% -highest overlap 89%) ( Figure 3B).
Next, we compared identified miRNA gene targets in islets with genes that showed altered expression in T2D donor islets in published bulk-and single-cell RNA-seq data (Bugliani et al., 2013;Fadista et al., 2014;Gunton et al., 2005;Lawlor et al., 2017;Segerstolpe et al., 2016;Solimena et al., 2018;Taneera et al., 2012;Xin et al., 2016). We decided to include all 3,365 unique genes with altered expression in T2D from all studies (Table S7), as there was only a small overlap between the differentially expressed genes across the studies. This could be due to differences in RNA-profiling technologies (bulk sequencing, single-cell sequencing, and microarrays) and parameters (e.g. islet isolation procedure, sample number, statistical analysis, and variable control) of each study (Table S7). The overlap between all validated targets genes of the two datasets, as well as the proportion of these that are islet specific and differently expressed in T2D, is shown in Figure 3B. It is worth mentioning that iScience Article miR-130b-3p in the Up-Pos dataset and miR-23b-3p in the Down-Neg dataset had the highest number of validated and expressed gene targets in islets. The targets of these two miRNAs also displayed the highest overlap with previously reported T2D differently expressed genes.
miRNA gene target clustering and functional annotation reveals two islet-specific clusters Next, we questioned whether miRNAs in the Up-Pos and Down-Neg datasets are implicated into processes related to b-cell function. We first clustered the genes into high-correlating groups (clusters) by performing a weighted gene correlation network analysis (WGCNA), as we wanted to identify gene co-expression networks regulated by miRNAs. This allowed us to divide the genes into 39 and 43 clusters of various sizes in the Up-Pos and Down-Neg datasets, respectively (Table S10). After summarizing the expression profile of the individual clusters using their eigengenes, which represent the first principal component, we were able to correlate each cluster with phenotypic traits such as HbA1c, diabetic status, stimulatory index, human islet purity, and the normalized expression (log2) of the insulin and glucagon genes. The selection of the expression levels of insulin and glucagon as traits enabled us to determine how a-or b-cell-enriched the clusters were. The Up-Pos dataset had 11 b-cell-enriched clusters, six a-cell-enriched clusters, and 17 clusters enriched for both a-and b-cells, while the Down-Pos dataset had 14 b-cell-enriched clusters, 14 a-cell-enriched clusters, and 13 enriched for both. As explained below, the two islet specific clusters which were further investigated were enriched for both a-and b-cell phenotypic traits ( Figure 4A).
Enrichment analysis of functional annotation terms was then performed to associate the different clusters to biological pathways. The miRNA gene targets of each cluster were scanned for enriched Gene Ontology (GO) terms and terms belonging to the Reactome and KEGG pathway databases (Tables S8-S9). One cluster in the Up-Pos dataset (Up-Pos_Cluster 3) and one cluster in the Down-Neg dataset (Down-Neg_Cluster 1) included enriched terms/pathways related to a-and b-cells. In the Up-Pos dataset, we observed terms relative to exocytosis (''exocytotic vesicle'', ''vesicle docking'') and ''response to insulin'', while in the Down-Neg dataset we came across the term ''insulin secretion'' in all three databases, as well as pathways such as ''insulin processing'', ''regulation gene expression'', and ''development in b-cells'', and terms related to responses to insulin and glucagon signaling pathways ( Figure 4B).
Next, we compared all 82 clusters with the list of previously reported genes with altered expression in T2D (Table S7). Interestingly, Up-Pos_Cluster 3 and Down-Neg_Cluster 1 included a higher proportion of islet T2D genes (Up-Pos_Cluster 3: 42%, Down-Neg_Cluster 1: 40%) compared to other clusters of similar size in each dataset ( Figure 4C, Table S10). Finally, we observed a broad specificity of differentially expressed miRNAs ( Figure 2) regarding gene targeting in the different clusters, as almost all miRNAs in each dataset target genes belonging to medium-sized and large clusters ( Figure S3).

Genetic regulation of altered and abundant miRNAs
eQTLs can provide an important step to link the miRNAs with T2D-risk variations derived from previous GWAS studies Piñ ero et al., 2020;Ramos et al., 2014). We therefore sought to determine cis-eQTLs mapping to within a 1Mb window of the starting position of each miRNA for our miRNA sets. This analysis was performed in three different miRNA datasets: the Up-Pos, Down-Neg, and the 15 most abundant miRNAs (top15). No eQTL was found to be statistically significant after p value correction for multiple testing which could be explained by the low study power of the analysis (n=18). Given that multiple nominal signals at the same locus are suggestive of a link between the genetic variant and miRNA expression, we considered eQTLs showing a nominal p value<0.05 (Table S11).
Co-localization of the miRNA eQTLs with T2D-risk loci would provide additional hints of the implications of these miRNAs in the disease. After assembling GWAS data from multiple sources, we generated a list of ll OPEN ACCESS 6 iScience 25, 103995, April 15, 2022 iScience Article SNPs that have been associated with T2D and related glycemic and insulin secretion traits (Methods and Table S12). The list was used to check for overlaps between these SNPs and the miRNA eQTLs in the three miRNA datasets. In total, three miRNAs (miR-7-5p, miR-126-3p, miR-1236) had nominal eQTLs overlapping with suggestive signals for T2D, two miRNAs (miR-130b-5p, miR-1275) had nominal eQTLs overlapping with suggestive signals for insulin secretion indices, and miR-194-5p had nominal eQTLs overlapping with both signals ( Figure 5, Table S13).

DISCUSSION
In this study, we investigated miRNAs-mRNA networks in pancreatic islets and their potential role in the pathogenesis of T2D. We identified miRNAs with altered expression in hyperglycemia, as well as their target genes. We also explored clusters of islet miRNA-mRNA networks and identified specific clusters associated with insulin and glucagon expression as well as T2D. Moreover, we found suggestive eQTLs iScience Article linking genetic variation to the differentially expressed miRNAs, which overlapped with T2D and insulin secretion GWAS signals.
The profiling of the miRNA expression levels in human islets revealed a limited number of highly expressed miRNAs, some of which have been previously described to be important for islet function. Among the 15 most expressed miRNAs, we identified miR-375, miR-29 family (miR-29a/b/c), members of miR-200 family (miR-200c/miR-141), and miR-7. miR-375 is important in the regulation of many functions in the b-cell including proliferation, insulin synthesis, processing, and secretion (Eliasson, 2017). Members of the miR-29 and miR-200 families are implicated in b-cell apoptosis and insulin secretion, and miR-7 is involved in insulin secretion and development (reviewed in (Eliasson and Esguerra, 2020;LaPierre and Stoffel, 2017)). The most abundant miRNA in our assay, miR-1246, has not previously been described in a b-cell context. However, it has been classified as a potential serum biomarker for diabetes (Vasu et al., 2019) and pancreatic cancer prognosis (Wei et al., 2020).
We identified 37 upregulated and 26 downregulated miRNAs in islets from IGT/T2D donors compared to NGT. Other studies have previously shown a similar pattern regarding the larger number of upregulated miRNAs in diabetes/obesity in human, mice, and rat models Locke et al., 2014;Nesca et al., 2013;Zhao et al., 2009). Hemoglobin A1c (HbA1c) is used as a measure of the average blood glucose levels over the last 2-3 months. As a means to increase our confidence that the expression of the differentially expressed miRNAs is associated with the donor glycemic state, we selected those that also showed a direct correlation with HbA1c in the appropriate direction. Thus, we ended up with a set of nine miRNAs upregulated and positively correlated with HbA1c (Up-Pos) and a set of 22 miRNAs downregulated and negatively correlated with HbA1c (Down-Neg). Interestingly, only five out of the 31 differentially expressed miRNAs are among the most abundant ones (miR-141-3p, miR-21-5p, miR-29b-3p, miR-24-3p, and miR-7- Figure 5. eQTL locus zoom plots of miRNAs with suggestive signals of association with insulin secretion/T2D-risk loci The points represent eQTLs within 1Mb of the start of the corresponding miRNA and are plotted against the significance of their correlation to the expression of the miRNA (-log10 p value). eQTLs above the red line have a nominal p value<0.05. Highlighted are the eQTLs that overlap T2D-risk variants or variants associated with insulin secretion traits. eQTLs in very close proximity may overlap, so it is recommended to also consider the complete list in Table  S13. iScience Article 5p; all downregulated), meaning that regulation and abundance do not necessarily coincide. Moreover, we found miRNAs in the Up-Pos and Down-Neg datasets that have been investigated in separate studies with a defective b-cell function perspective (Eliasson and Esguerra, 2020;LaPierre and Stoffel, 2017).
Instead of studying the effects of individual miRNAs on the b-cell functionality, we implemented a network approach, in which the whole list of differentially expressed miRNAs could point to potentially dysregulated pathways in donors with altered glycemic status. The first step was to match the miRNAs of both sets with validated gene targets from public sources. We focused on validated rather than predicted miRNA targets as, despite considerable advances in the field of miRNA target prediction, the predictive algorithms are still far from perfect due to the inconsistency in miRNAs pairing with their target genes (Helwak et al., 2013;Liu and Wang, 2019). The fact that the list of the miRNA with most targets (miR-23b-3p) includes 3,529 genes is a good example of how a few thousand miRNAs can regulate over 60% of protein-coding genes (Friedman et al., 2009).
After filtering out gene targets not expressed in human islets, we implemented WGCNA to identify clusters of highly correlated miRNA gene targets. Functional annotation revealed two interesting clusters, one in the Up-Pos and one in the Down-Neg dataset, that most likely constitute islet-specific networks as they were positively associated with islet purity, insulin, and glucagon expression levels. This was also supported by the enrichment of functional terms related to islet function, such as ''insulin biosynthesis and secretion'', ''b-cell development'', and ''response to insulin''. Of interest was Down-Neg_Cluster 1, comprising 1,019 correlated genes in total. Among these, there were genes crucial for the b-cell development and function, such as NEUROD1 (Gu et al., 2010), PAX6 (Gosmain et al., 2012), NKX2-2 (Doyle and Sussel, 2007), FOXO1 (Buteau and Accili, 2007), and KCNB1 (Fu et al., 2017). Hence, gene interactions can produce a highly complex network that fine-tunes distinct and vital functions within the islet. However, only some of the genes in the Down-Neg_Cluster 1 was previously reported to be differentially expressed in T2D (410/1,019 genes) (Bugliani et al., 2013;Fadista et al., 2014;Gunton et al., 2005;Lawlor et al., 2017;Segerstolpe et al., 2016;Solimena et al., 2018;Taneera et al., 2012;Xin et al., 2016) This can be explained by the disease module hypothesis, wherein a functional gene association will not fundamentally correspond to a specific disease phenotype (Barabá si et al., 2011), meaning that dysregulation of only some of the genes in the module (cluster) can potentially be implicated in the pathophysiology of T2D.
GWAS facilitate the identification of disease-risk loci through associating genetic variants to diseases such as T2D. These loci are named after the nearest gene, which may not necessarily be the causal gene. In order to identify the causal gene, these genetic variants are investigated for their association with gene or miRNA expression, so called expression quantitative trait loci (eQTLs). Mapping miRNA expression to proximal loci did not result in any significant associations, potentially due to the small power of our study (n=18). However, even larger-scale studies did not manage to identify statistically significant eQTLs. In a previous study with 176 lymphoblastoid cell lines, no cis-miR-eQTLs were found (Gamazon et al., 2012), while another study in primary fibroblasts of 180 new-borns revealed only 12 cis-miR-eQTLs with FDR<0.5 (Borel et al., 2011). The complications of correlating SNPs to miRNA is confirmed by a large-scale study, which reports that single cis-miRNA eQTLs can explain only a small portion of the expression variability of their associated miRNAs (1.3%) compared to the cis-mRNA eQTLs, which explain 33%-53% of variance in the expression levels of their associated mRNAs (Huan et al., 2015). This implies that there might be more suitable statistical models, other than linear correlations, that could identify cis-miRNA-eQTLs with more precision.
However, even multiple suggestive signals linking genetic variation to miRNA expression, which in turn is correlated to a relevant phenotype, can be indicative of their involvement in that particular phenotype. The overlap between 13 miRNA eQTLs with T2D loci and insulin secretion indices suggests the implication of differentially expressed miRNA regulating genes of Up-Pos and Down-Neg datasets in T2D. Specifically, eQTL of miR-194 overlaps with suggestive signals for both T2D and insulin secretion phenotypes. This miRNA has previously been shown to be a biomarker for diabetes incidence (Jaeger et al., 2018), as well as a regulator of glucose metabolism in skeletal muscle (Latouche et al., 2016). Another interesting miRNA, miR-126, whose eQTLs showed suggestive signals for T2D, was found to be a potent biomarker for early prognosis of diabetes (Liu et al., 2014)  iScience Article diabetes'' (Piñ ero et al., 2020), an interesting finding considering the role of miR-1236 as a biomarker for latent autoimmune diabetes in adults (LADA) (Yu et al., 2019). Out of the two miRNAs whose eQTL profile suggests signals for insulin secretion traits, miR-130b has been shown to affect intracellular ATP levels in the pancreatic b-cell, indicating a possible effect on insulin secretion (Ofori et al., 2017), while miR-1275 has not been associated with T2D or insulin secretion before. Notably, none of the 13 eQTLs with suggestive T2D-risk/insulin secretion signals has been functionally linked to T2D so far, making them suitable candidates for future T2D studies. Moreover, additional eQTLs of the six miRNAs with strong nominal signals that have not been associated with any existing traits before could be considered potential targets of interest in the pursuit of a more complete comprehension of T2D etiology.
In conclusion, our work demonstrates the complexity of the miRNA-mRNA network regulation in pancreatic islets. Therefore, efforts toward miRNA-based therapeutic strategies should focus on the comprehension of the role of individual miRNAs as network components rather than individual gene modulators.

Limitations of the study
The differentially regulated miRNAs identified in this study do not overlap with those described in other studies that compare healthy and T2D human islets (Kameswaran et al., 2014;Locke et al., 2014). The variability of individual donors in terms of their underlying pathophysiological condition, the differences in islet isolation and culturing techniques, the small number of samples under comparison and the distinct profiling methods could all provide plausible explanations for this discrepancy. Another caveat is the small number of samples included in the study, a common problem due to the limited supply of human islets. This is why the principal aim of this study is to explore the complexity of the miRNA-mRNA networks and how they can potentially contribute to T2D onset, rather than defining casual miRNA-mRNA-phenotype relationships. Moreover, miRNA-regulated network analysis in islets was based on the collection of experimentally validated miRNA gene target data from multiple studies. However, it should be noted that despite revealing genuine miRNA-mRNA interactions, high-throughput miRNA-target validation techniques are not necessarily accompanied by functional target regulation (Liu and Wang, 2019) and results should be interpreted with caution.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

ACKNOWLEDGMENTS
We gratefully acknowledge Leif Groop for establishing and mentoring the islet Geneview project from which the human islet data are sourced. We thank Anna-Maria Veljanovska Ramsay at the Department of Clinical Sciences Malmö , Lund University for technical assistance. This work was supported by a grant from the Swedish Foundation for Strategic Research, Sweden (IRC-LUDC; DnrIRC15-0067) and the ll OPEN ACCESS